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We present a comprehensive theoretical analysis and computational study of four-wave mixing 
(FWM) of optical pulses co-propagating in one-dimensional silicon photonic crystal waveguides (Si- 
PhCWGs). Our theoretical analysis describes a very general set-up of the interacting optical pulses, 
namely we consider nondegenerate FWM in a configuration in which at each frequency there exists 
a superposition of guiding modes. We incorporate in our theoretical model all relevant linear optical 
effects, including waveguide loss, free-carrier (FC) dispersion and FC absorption, nonlinear optical 
effects such as self- and cross-phase modulation (SPM, XPM), two-photon absorption (TPA), and 
cross-absorption modulation (XAM), as well as the coupled dynamics of FCs and optical field. In 
particular, our theoretical analysis based on the coupled-mode theory provides rigorously derived 
formulae for linear dispersion coefficients of the guiding modes, linear coupling coefficients between 
these modes, as well as the nonlinear waveguide coefficients describing SPM, XPM, TPA, XAM, 
and FWM. In addition, our theoretical analysis and numerical simulations reveal key differences 
between the characteristics of FWM in the slow- and fast-light regimes, which could potentially 
have important implications to the design of ultra-compact active photonic devices. 

PACS numbers: 78.67.Pt, 78.20.Bh, 42.65.Wi, 42.70.Qs, 42.65.Ky 


I. INTRODUCTION 


One of the most promising applications of photonics 
is the development of ultra-compact optical intercon¬ 
nects for chip-to-chip and even intra-chip communica¬ 
tions. The driving forces behind research in this area 
are the perceived limitations at high frequency of cur¬ 
rently used copper interconnects [l|, combined with a 
rapidly increasing demand to move huge amounts of data 
within increasingly more confined yet increasingly intri¬ 
cate communication architectures. An approach show¬ 
ing great potential towards developing optical intercon¬ 
nects at chip scale is based on high-index contrast opti¬ 
cal waveguides, such as silicon photonic waveguides (Si- 
PhWGs) implemented on the silicon-on-insulator mate¬ 
rial platform ii- Among key advantages provided by 
this platform are the increased potential for device in¬ 
tegration facilitated by the enhanced confinement of the 
optical field achievable in high-index contrast photonic 
structures, as well as the particularly large optical nonlin¬ 
earity of silicon, which makes it an ideal material for ac¬ 
tive photonic devices. Many of the basic device function¬ 
alities required in networks-on-chip have in fact already 
been demonstrated using Si-PhWGs, includ ing pa ramet- 
ric amplification J0-Q , optical modulation (91411 , p ulse 
compression [l^, [l 3 1 supercontinuum generation 3 3 ; 
pulse self-steepening [l6|, modulational instability 17|, 
and four-wave mixing (FWM) [3423; for a review of 
optical properties of Si-PhWGs see [3l- However, since 
the parameter space of Si-PhWGs is rather limited, there 
is little room to engineer their optical properties. 

A promising solution to this problem has its roots in 
the advent of photonic crystals (PhCs) in the late 80’s 
Thus, by patterning an optical medium in a 


periodic manner, with the spatial periods of the pattern 
being comparable to the operating optical wavelength, 
the optical properties of the resulting medium can be 
modified and engineered to a remarkable extent. Fol¬ 
lowing this approach, a series of photonic devices have 
been demonstrated u sing PhGs, including optical waveg¬ 
uides and bends 126143011 . optical micro-cavities |3ll435|| . 
and optical filters (3; 1^ • One of the most effective ap¬ 
proaches to affect the optical properties of PhCs is to 
modify the group-velocity (GV), Vg, of the propagating 
modes. Unlike the case of waves propagating in regular 
optical media, whose GV can hardly be altered, by vary¬ 
ing the geometrical parameters of PhCs one can tune the 
corresponding GV over many orders of magnitude. Per¬ 
haps the most noteworthy implication of the existence 
of optical modes with significantly reduced GV, the so- 
called slow-light [ssl - ldoll , is that both linear and nonlinear 
optical effects can be dramatically enhanced in the slow- 
light regime [4ll - l4l| . 

One of the most important nonlinear optical process, 
as far as nonlinear optics applications are concerned, is 
FWM. In the generic case, it consists of the combina¬ 
tion of two photons with frequencies, wi and W 2 , belong¬ 
ing to two pump continuous-waves (CWs) or pulses, fol¬ 
lowed by the generation of a pair of photons with fre¬ 
quencies UJ 3 and W 4 . The energy conservation requires 
that LOi ijj 2 = LOz In practice, however, an easier 

to implement FWM configuration is usually employed, 
namely degenerate FWM. In this case one uses just one 
pump with frequency, Wp, the generated photons belong¬ 
ing to a signal {oJs) and an idler (oji) beam; in this case 
the conservation of the optical energy is expressed as: 
2ujp = ujg+uji. Among the most important applications of 
degenerate FWM it is noteworthy to mention optical am- 
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FIG. 1. (a) Geometry of the ID Si-PhC slab waveguide. The 
height of the slab is h = 0.6a and the radius of the holes is 
r = 0.22a. The primed coordinate system shows the principal 
axes of the Si crystal with the input facet of the waveguide 
in the (110) plane of the Si crystal lattice, (b) Projected 
band structure. Dark yellow and brown areas correspond to 
slab leaky and guiding modes, respectively. The red and blue 
curves represent the guiding modes of the ID waveguides. 


plification, wavelength generation and conversion, phase 
conjugation, generation of squeezed states, and supercon¬ 
tinuum generation. While FWM has been investigated 
theoretically and experimentally in PhC waveguides [4§- 
[ 5 ^ and long-period Bragg waveguides [ 2 ^ [Hj, a com¬ 
prehensive theory of FWM in silicon PhC waveguides 
(Si-PhCWGs), which rigorously incorporates in a uni¬ 
tary way all relevant linear and nonlinear optical effects 
as well as the influence of photogenerated free-carriers 
(FCs) on the pulse dynamics is not available yet. 

In this article we introduce a rigorous theoretical model 
that describes FWM in Si-PhCWGs. Our model cap¬ 
tures the influence on the FWM process of linear optical 
effects, including waveguide loss, FC dispersion (FCD) 
and FC absorption (FCA), nonlinear optical effects such 
as self- and cross-phase modulation (SPM, XPM), two- 
photon absorption (TPA), and cross-absorption modula¬ 
tion (XAM), as well as the mutual interaction between 
FCs and optical field. We also illustrate how our model 
can be applied to investigate the characteristics of FWM 
in the slow- and fast-light regimes, showing among other 
things that by incorporating the effects of FCs on the 
optical pulse dynamics new physics emerge. One note¬ 
worthy example in this context is that the well-known 
linear dependence of FCA on v~^ is replaced in the slow- 
light regime by a v~^ power-law dependence. 

The remaining of the paper is organized as follows. In 
the next section we present the optical properties of the 
PhC waveguide considered in this work. Then, in Sec. Ell 
we develop the theory of pulsed FWM in Si-PhCWGs 
whereas the particular case of degenerate FWM is ana¬ 
lyzed in Sec. |lVl Then, in Sec. El we apply these theo¬ 


retical tools to explore the physical conditions in which 
efficient FWM can be achieved. The results are subse¬ 
quently used, in Sec. ED to study via numerical simula¬ 
tions the main properties of pulsed FWM in Si-PhCWGs. 
We conclude our paper by summarizing in the last section 
the main findings of our article and discussing some of 
their implications to future developments in this research 
area. Finally, an averaged model that can be used in the 
case of broad optical pulses is presented in an Appendix. 


II. DESCRIPTION OF THE PHOTONIC 
CRYSTAL WAVEGUIDE 

In this section we present the geometrical and material 
properties of the PhC waveguide considered in this work, 
as well as the physical properties of its optical modes. 
Thus, our Si-PhCWG consists of a one-dimensional (ID) 
waveguide formed by introducing a line defect in a two- 
dimensional (2D) honeycomb-type periodic lattice of air 
holes in a homogeneous slab made of silicon (a so-called 
W1 PhC waveguide). The line defect is oriented along 
the z-axis, which is chosen to coincide with one of the 
FAT symmetry axes of the crystal, and is created by fill¬ 
ing in a row of holes [see Fig. [I](a)]. The slab height is 
h = 0.6a and the radius of the holes is r = 0.22a, where 
a = 412 nm is the lattice constant, whereas the index of 
refraction of silicon is nsi = n = 3.48. 



FIG. 2. Left (right) panels show the amplitude of the normal¬ 
ized magnetic field of the y-odd (y-even) mode, calculated 
in the plane a; = 0 for five different values of the propagation 
constant, kz- From top to bottom, the panels correspond to 
the Bloch modes indicated in Fig.[T](b) by the circles A, B, C, 
D, and E, respectively. 
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FIG. 3. (a), (b), (c), and (d) Frequency dependence of waveg¬ 
uide dispersion coefficients Ug, 132, Pi, and Pi, respectively, de¬ 
termined for the even and odd modes. Light green, blue, and 
red shaded regions correspond to slow-light regime, defined 
as Ug > 20. The dashed vertical line in panel (b) indicates 
the zero-GVD wavelength. 


The defect line breaks the discrete translational sym¬ 
metry of the photonic system along the y-axis, so that the 
optical modes of the waveguide are invariant only to dis¬ 
crete translation along the 2 -axis [s^. Moreover, based 
on experimental considerations, we restrict our analysis 
to in-plane wave propagation, namely the wave vector, 
k, lies in the a; = 0 plane. The component, on the 
other hand, can be restricted to the first Brillouin zone, 
kz S [—7r/a, 7r/a], which is an immediate consequence of 
the Bloch theorem. Under these circumstances, we de¬ 
termined numerically the photonic band structure of the 
system and the guiding optical modes of the waveguide 
using MPB, a freely available code based on the plane- 
wave expansion (PWE) method [s^. To be more specific, 
we used a supercell with size of 6a x 19-\/3/2a x a along 
the X-, y-, and 2-axis, respectively, the corresponding step 
size of the computational grid being a/60, a-\/3/120, and 
a/60, respectively. Figure [TJb) summarizes the results 
of these calculations. Thus, the waveguide has two fun¬ 
damental TE-like optical guiding modes located in the 
band-gap of the unperturbed PhC, one y-even and the 
other one y-odd. 

In order to better understand the physical properties of 
the optical guiding modes, we plot in Fig.[2]the profile of 


the magnetic field H^, which is its only nonzero compo¬ 
nent in the a; = 0 symmetry plane. These field profiles, 
calculated for several values of k^, show that although 
the optical field is primarily confined at the location of 
the defect (waveguide), for some values of kz it is rather 
delocalized in the transverse direction. This field delo¬ 
calization effect is particularly strong in the spectral do¬ 
mains where the modal dispersion curves are relatively 
flat, namely in the so-called slow-light regime, and in¬ 
creases when the group index of the mode, defined as 
Ug = c/vg, increases. 

The dispersion effects upon pulse propagation in the 
waveguide are characterized by the waveguide dispersion 
coefficients, defined as /3„ = <F‘kz/duj'^. In particular, the 
first-order dispersion coefficient is related to the pulse GV 
via Pi = 3/vg, whereas the second-order dispersion coef¬ 
ficient, P 2 , quantifies the GV dispersion (GVD) as well as 
pulse broadening effects. The wavelength dependence of 
the first four dispersion coefficients, determined for both 
guided modes, is presented in Fig.O the shaded areas in¬ 
dicating the spectral regions of slow-light. For the sake of 
clarity, we set the corresponding threshold to c/vg = 20, 
that is the slow-light regime is defined by Ug > 20. As it 
can be seen in Fig. |3l the even mode possesses two slow- 
light regions, one located at the band-edge (A ~ 1.6 pm) 
and the other one at k^ ~ 0.3(27r/a), i.e. A ~ 1.52pm, 
whereas the odd mode contains only one such spectral 
domain located at the band-edge (A ~ 1.67 pm). More¬ 
over, the even mode can have both positive and negative 
GVD, the zero-GVD point being at A = 1.56 pm, whereas 
the odd mode has normal GVD {P 2 > 0) throughout. 
Since usually efficient FWM can only be achieved in the 
anomalous GVD regime {P 2 < 0), we will assume that 
the interacting pulses propagate in the even mode unless 
otherwise is specified. 


III. DERIVATION OF THE MATHEMATICAL 
MODEL 

This section is devoted to the derivation of a system 
of coupled-mode equations describing the co-propagation 
of a set of mutually interacting optical pulses in a Si- 
PhGWG, as well the influence of photogenerated FCs 
on the pulse evolution. We will derive these coupled¬ 
mode equations in the most general setting, namely the 
nondegenerate FWM, then show how they can be applied 
to a particular case most used in practice, the so-called 
degenerate FWM configuration. Our derivation follows 
the general approach used to develop a theoretical model 
for pulse propagation in silicon waveguides with uniform 
cross-section and Si-PhGWGs [^ . 


A. Optical modes of photonic crystal waveguides 

In the presence of an external perturbation described 
by the polarization, Ppert(r,a;), the electromagnetic field 
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of guiding modes with frequency, w, is described by the 
Maxwell equations, which in the frequency domain can 
be written in the following form: 

V X E(r, w) = za;^H(r,w), (la) 

V X H(r,w) =-iw[ec(r,w)E(r,u;)+Ppert(r,w)], (lb) 

where ^ is the magnetic permeability, which in the case 
of silicon and other nonmagnetic materials can be set to 

= fio, ec(r, w) is the dielectric constant of the PhC, and 
E and H are the electric and magnetic fields, respectively. 
In our case, Ppert is the sum of polarizations describing 
the refraction index change induced by photogenerated 
FCs and nonlinear (Kerr) effects. 

In order to understand how the modes of the PhC 
waveguide are affected by external perturbations, let us 
consider first the unperturbed system, that is Ppert = 0. 
Thus, let us assume that, at the frequency w, the un¬ 
perturbed PhC waveguide has M guiding modes. It fol¬ 
lows then from the Bloch theorem that the fields of these 
modes can be written as: 

Em,<r(r,w) = e„j^(r,w)e*'^^'""', m = 1, 2 ,..., M, (2a) 

Hr„,(r,u;) =h„<r(r,cc)e*"^"*", m = l,2,...,M, (2b) 

where Pm is the mth mode propagation constant and 
cr = -I- (cr = —) denotes forward (backward) propagating 
modes. Here, we consider that the harmonic time depen¬ 
dence of the fields was chosen as The mode am¬ 

plitudes Bma and hma are periodic along the z-axis, with 
period a. Moreover, the forward and backward propa¬ 
gating modes obey the following symmetry relations: 

e^_(r,a;) = e);^_p(r,w), (3a) 

h,„_(r,w) = -h);,^(r,a;), (3b) 

where the symbol denotes complex conjugation. As 
such, one only has to determine either the forward or the 
backward propagating modes. 

The guiding modes can be orthogonalized, the most 
commonly used normalization convention being 


are the electric and magnetic energy of the mode, respec¬ 
tively, and 14eii is the volume of the unit cell. Note that 
in Eq. ([5]) we used the fact that the mode contains equal 
amounts of electric and magnetic energy. 

It should be stressed that the waveguide modes defined 
by Eqs. ([2|) are exact solutions of the Maxwell equations 
(HD with Ppert = 0 , and thus they should not be confused 
with the so-called local modes of the waveguide. The lat¬ 
ter modes correspond to waveguides whose optical prop¬ 
erties vary adiabatically with z, on a scale comparable to 
the wavelength and have been used to describe, e.g., wave 
propagation in tapered waveguides or p ulse propa¬ 
gation in ID long-period Bragg gratings [^. 


B. Perturbations of the photonic crystal waveguide 


Due to the photogeneration of FCs and nonlinear opti¬ 
cal effects, the dielectric constant of Si-PhCWGs under¬ 
goes a certain local variation, (5e(r), upon the propagation 
of optical pulses in the waveguide. The corresponding 
perturbation polarization, Ppert in Eq. (Ilbl) . can be di¬ 
vided in two components according to the physical effects 
they describe: the linear change of the dielectric constant 
via generation of FCs and the nonlinearly induced vari¬ 
ation of the index of refraction. 

Assuming an instantaneous response of the medium, 
the linear contribution to Ppert, <5Piin(r,t), is written as: 


(5Piin(r, t) = [detc(r) + (5eioss(r)] E(r, t), 


where 


detc(r) = (2eonSn{c 

^eioss(r) = ^-QiinS(r). 

UJ 


.eocn 
• I - Otic 


) S(r), 


(7) 

(8a) 

(8b) 


Here, is the intrinsic loss coefficient of the waveg¬ 
uide and E(r) is the characteristic function of the do¬ 
main where FCs can be generated, namely E = 1 in the 
domain occupied by Si and E = 0 otherwise. Based on 
the Drude model, the FC-induced change of the index of 
refraction, dntc, and FC losses, ate, are given by M- 


I (^^7n<7 ^ ^m'er' ^ 

J s 

(4) 

where Pm is the power carried by the mth mode. This 
mode power is related to the mode energy contained in 
one unit cell of the PhC waveguide, Wm, via the relation: 

wel wmag 2W‘^^ 

^ + ^ (5) 

a a a 

where 

" 7 / ^(a;ec)|e™„(r,a;)pdH, (6a) 

WC® = T / Mo|h„^(r,w)pdH, (6b) 


(9a) 
(9b) 

Here, e is the charge of the electron, /ie (/x/j) is the elec¬ 
tron (hole) mobility, nice = 0.26mo {nich = 0.39mo) is 
the conductivity effective mass of the electrons (holes), 
with mo the mass of the electron, and Ne (Nh) is the in¬ 
duced variation of the electrons (holes) density (in what 
follows, we assume that Ne = Nh = N). 

The nonlinear contribution to Ppert, <^Pni(r,t), is de¬ 
scribed by a third-order nonlinear susceptibility, y;*^^)(r), 
and can be written as: 

dPni(r,t) = eoX^^^(r):E(r,t)E(r,t)E(r,t). (10) 
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The real part of the susceptibility describes para¬ 
metric optical processes such as SPM, XPM, and FWM, 
while the imaginary part of corresponds to TPA and 
XAM. Note that in this study we neglect the stimulated 
Raman scattering effect as it is assumed that the frequen¬ 
cies of the interacting pulses do not satisfy the condition 
required for an efficient, resonant Raman interaction. 

Since silicon belongs to the crystallographic point 
group m3m the susceptibility tensor has 21 nonzero 
elements, of which only 4 are independent, namely, Xim, 
X1122, X1212, and X1221 [i^- In addition, the frequency 
dispersion of the nonlinear susceptibility can be neglected 
as we consider optical pulses with duration of just a few 
picoseconds or larger. As a consequence, the Kleinman 
symmetry relations imply that X1122 = X1212 = Xi22i- 

-s ( 3 ) 

Moreover, experimental studies have shown that Xiiii = 
2.36 xJ^ 22 within a broad frequency range. There¬ 
fore, the nonlinear optical effects considered here can be 
described by only one element of the tensor x^^^- 

Because of fabrication considerations, in many in¬ 
stances the waveguide is not aligned with any of the 
crystal principal axes and as such these axes are different 
from the coordinate axes in which the optical modes are 
calculated. Therefore, one has to transform the tensor 
from the crystal principal axes into the coordinate 
system in which the optical modes are calculated (s^ . 

Xijkl = RiaRjpRkjRlSXap^s, ( 11 ) 

where x^^^^ is the nonlinear susceptibility in the crystal 
principal axes and R is the rotation matrix that trans¬ 
forms one coordinate system into the other. In our case, 
R is the matrix describing a rotation with 7 r /4 around 
the a;-axis (see Fig.[T|). 


Let us consider now a nondegenerate FWM process in 
which two pulses at carrier frequencies Qi and (D 2 interact 
and generate two optical pulses at carrier frequencies (Da 
and W 4 , with the energy conservation expressed as wi -I- 
W 2 = W 3 -|- W 4 . Then, in the Lorentz reciprocity theorem 
given by Eq. CSl) we choose as the first set of fields a 
mode of the unperturbed waveguide (Ppert = 0), which 
corresponds to the frequency to a = oJi, where uJi is one of 
the carrier frequencies wi, W 2 , W 3 , or 0)4: 


Ea{r,uji) 


Ha(r,a;*) 


ip j 

hniPi(r,Wi) ip. 


(14a) 


(14b) 


where pi = ±1 and Ui is an integer, 1 < < W, * = 

1,..., 4, with Ni being the number of guiding modes at 
the frequency Wi. In Eqs. d, and in what follows, a bar 
over a symbol means that the corresponding quantity is 
evaluated at one of the carrier frequencies. 

As the second set of fields we take those that propagate 
in the perturbed waveguide, at the frequency oJb = w. 
These fields are written as a series expansion of the guid¬ 
ing modes at frequencies w,, z = 1,..., 4, thus neglecting 
the frequency dispersion of the guiding modes and the ra¬ 
diative modes that might exist at the frequency w. This 
approximation is valid as long as all interacting optical 
pulses have narrow spectra centered at the corresponding 
carrier frequencies, that is the physical situation consid¬ 
ered in this work. In particular, this modal expansion 
becomes less accurate when any of the pulses propagates 
in the slow-light regime, as generally the smaller the GV 
of a mode is the larger its frequency dispersion is. Thus, 
the second set of fields are expanded as: 


C. Coupled-mode equations for the optical field 


In order to derive the system of coupled-mode equa¬ 
tions describing pulsed FWM in Si-PhCWGs we employ 
form of the Lorentz reciprocity theorem 
) this end, let us consider two solutions 
of the Maxwell equations (U), [Ea(r, Wa), Ha(r, Wa)] and 
[Eb(r, Wb),H{,(r,a;b)], which correspond to two different 
spatial distribution of the dielectric constant, ea(r,uJa) 
and e{,(r, Wb), respectively. If we insert the vector E, de¬ 
fined as E = E{, X H* -I- E* X Hf,, in the integral identity: 


the co njug ated 

[13, [M^. T 


V •FdS = — 

' oz, 


F-MS- 


E • ndL 


dS 


( 12 ) 


where S is the transverse section at position, z, and dS 
is the boundary of 5, and use the Maxwell equations, we 
arrive at the following relation: 


dz 


UbdS 


f F-MS = tpoiLOb-u;a) f H: 

Js Js 

+ i / (wbCb - uJaea)Fl ■ FbdS - f F ■ ndl. 
JS JdS 


( 13 ) 


Fb{r,u;)=^Y. 

j=l rrijaj JPruij 

(15a) 

H ( \ \ ' \ ' (7) ( \ (^1 ^i) iajBm-Z 

= 2 _^ 2 ^ a)l>.^.{z,u) - ^-j ==—e . 

j=lmjaj y Rrrij 

(15b) 

With the fields normalization used in Eqs. m, the mode 
amplitudes amia-i{z,uj), i = 1,...,4, are measured in 
units of a/W. Note that since the optical pulses are 
assumed to be spectrally narrow, the mode amplitudes 
OmiCTi ( 2 :, w) have negligible values except when the fre¬ 
quency uj lies in a narrow spectral domain centered at 
the carrier frequency, uii. 

The dielectric constant in the two cases is Ca = ec(r, Coi) 
and £{, = ec(r, w)-l-(5e(r,a;), where ec{r,oj) is the dielectric 
constant of the unperturbed PhC. If the material disper¬ 
sion is neglected, ec(r,a;) = ec{r,uji) = ec(r). Inserting 
the fields given by Eqs. (03 and Eqs. (03 in Eq. (03 . and 
neglecting the line integral in Eq. (USD, which cancels for 
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exponentially decaying guiding modes, one obtains the 
following set of coupled equations: 


P^- 




dz 

(^) 

dz 


- 1 mj(7j 

^ i 


= K,pAU^)+ E ' 


luje 


4/^ 


jrrijc 


-'riipi 


Ppert ^ — 4;---;4, 

(16) 


where 


C%p,,m,a,{z) = 


i^rjjPm.-piP-nd^ 


Prii Pm j 

I (SmjCTj X hj^.p. + e„.p. X hmjtTj) ■ ids, 

Js 

(17a) 


TDl _ __ 

- 4p 


niPi,mjG j 


[^xo{uJ - a;i)|h„,p,p 

Ui Js 

+{ujec - ecWi)|en,p,P] dS, 

,: Pm J-Pi Prii)^ r 

{z) = - -/ [Po{u^ - w.) 


(17b) 




^^mjaj * ^riiPi P ^c^i'j^mjCTj ' ^riip,] dS. 

(17c) 

In Eq. (1161) and what follows a prime symbol to a sum 
means that the summation is taken over all modes, except 
that with j = i, mj = rii, and <7j = pi. Moreover, in 
deriving the Lh.s. of Eq. (ITCl) we used the orthogonality 
relation given by Eq. O- 

The time-dependent fields are obtained by integrating 
over all frequency components contained in the spectra 
of the system of interacting optical pulses: 

=) TE E 


j—l mj<7j 


Prr 


X + c.c. = \ [E(+)(r,t) -f E(-)(r,t) 


2 L 


TT(' ^ (?) ( \ (r, Wj) 


(18a) 


0 j — 1 TYljGj 


Prr 


X + c.c. = i [H(+)(r,t) + H(-)(r,t) 

(18b) 


where E(+)(r,t), H(+)(r,t) and E(“^(r,t), H^^)(r,t) are 
the positive and negative frequency parts of the spec¬ 
trum, respectively. 


Let us now introduce the envelopes of the interacting 
pulses in the time domain, An}pi{z,t), defined as the in¬ 
tegral of the mode amplitudes taken over the part of the 
spectrum that contains only positive frequencies, 

roo 

A^l^iz,t)= a%^{z,u:)e-^P-^^)*dw. (19) 

Jo 

With this definition, the time-dependent fields given in 
Eqs. become: 


E(r,t) =^E E 


j — 1 mjCTj 


H(r,t) = 2 E E 

j — 1 rrijCTj 


Prr 


+C.C.. (20b) 


Following the same approach, the time-dependent po¬ 
larization, too, can be decomposed in two components, 
which contain positive and negative frequencies, that is, 
it can be written as: 


Epert ( 1*7 ^) 



— p(+) 

ey pert 


pert(l',w)e “‘dw-fc.c. 

(r,0+P^(r70 • 


( 21 ) 


The next step of our derivation is to Fourier trans¬ 
form Eq. (HU) in the time domain. To this end, we first 
expand the coefficients Bh.p. and Dl^-p-^rnj< 7 j iii Taylor se¬ 
ries, around the carrier frequency tOi [note that according 
to Eq. (I17a|) . Clj.p. „^.^. is frequency independent]: 


B'n„. = E 


(Ai.,)" a»B; 


riipi 


g>i 


q\ duji 


E -o{<l)i 

q\ 

UJ—UJi 


nipi,mj 




(ACC,)^ 


g>i 


9 ! 


du}‘i 


■/o(9)y 


g>i 


9 ! 


(22a) 


(22b) 


where Awi = uj — Wi, i = 1,... ,4. Combining Eqs. p2a|) . 
(I17bl) . dS]), and dH) leads to the following expression for 
the dispersion coefficients, Pn}pi ■ 





SUz) 


& 


n —1 




g,n.i 


n>2, 


(23a) 

(23b) 
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where 




+ ^(t^ec)|e„,p,(r,Wi)|^ 


dS. 


(24) 


It can be easily seen from this equation that the aver¬ 
age of oil® lattice cell of the PhC waveguide 

is equal to 1, ie., 


1 pz+a 

a 


SUz')dz' = l. 


(25) 


Here and in what follows the tilde symbol indicates that 
the corresponding physical quantity has been averaged 
over a lattice cell of the waveguide. With this notation, 
Eqs. (1^ become: 


a;’;, ^«... = :J-. 


(26a) 


g,ni 


Qn-lni 

e. ^ n > 2. (26b) 


du} 


These relations show that /3i”p) = m ^^i® ii^h order 
dispersion coefficient of the waveguide mode character¬ 
ized by the parameters {ni,pi}, evaluated at w = Wi. 


We now multiply Eq. (HU) hy e and integrate 

over the positive-frequency domain. These simple calcu¬ 
lations lead to the time-domain coupled-mode equations 
for the field envelopes, An}pPz,t): 


Pi 


9^ 

dz 


T. Ecs 

J — 1 ^3 

3 7 ^ i 




Pi I'^^j ^j 


'dA% 

dz 


+i{crjPm,j PiPnpA^.^. — i 'y 


9>1 


q{q)i 

Ariipi 


dt 


A{i) 

■A^rn - 


riipi 


' Priipi, 


q'- 


*EE 

q'> 1 jrrij <7j 


Pi^nijCTj I ^0) 


dt 


•3 Pj 


i = 1, . . . ,4, 
(27) 


4\/ Pm 


’'riipi 


The temporal width of the pulses considered in this 
analysis is much smaller as compared to the nonlinear 
electronic response time of silicon and therefore the latter 
can be approximated to be instantaneous. In addition, 
we assume that the spectra of the interacting pulses are 
narrow and do not overlap. Under these circumstances, 
the optical pulses can be viewed as quasi-monochromatic 
waves and their nonlinear interactions can be treated in 
the adiabatic limit. Separating the nonlinear optical ef¬ 
fects contributing to the nonlinear polarization, one can 
express in the time domain this polarization as 67| : 


^ ^ 

5Pni,.D,(r,t) = ^ E 

mi(Ti 

0 ^ ^ 

+ 2 E E eoX^^nw*,-Wi,Wi):ep.e,(r,Wi)e;.g.(r,w,)e,„,a,(r,w.)l^p-e.P^ml 

Pi > rrii 




Prrii \/Prrii 

gidiPmiZ 

Ppi 'JPnni 


+ 9 E E ®0X^^^K,-Wj,Wi):ep^e,(r,a;j)e;.^^(r,a;j)e77».,7.(r,a;0|^^^^p74W^ 


j — 1 rriiGi 

3 # i Pj Bj 


oltyiPmA Z 


Ppj 'JPrrii 


+ 9 E 


Pjqkmi 

QjTk^l 


j[iBjPpj-'^kPq,,+<ZlPmi)z-{ujj-UJk+UJl)t] 


Ppi ^qk Pm I 


j ^ k ^ I ^ i 

Oj — Qk + Cji = uii 


(28) 


r 


This expression for the nonlinear polarization accounts represents SPM effects of the pulse envelopes, the second 
for the fact that the nonlinear susceptibility is invariant and third terms describe the XPM between modes with 
to frequency permutations. The first term in Eq. (1281) the same frequency and XPM between pulses propagat- 






















ing at different frequencies, respectively, whereas the last 
term describes FWM processes. 

If one inserts in Eq. (1271) the linear and nonlinear po¬ 
larizations given by Eq. © and Eq. (E51) . respectively, 
then discards the fast time-varying terms, one obtains 
the following system of coupled equations that governs 
the dynamics of the mode envelopes: 




^A{i) o{q)i 

Hriipi 


dz 


9>1 


q\ \ dt 


riipi 


E 


3 ( 9 )* 


q'- 




dt 




riipi ,miCT, 




3iuji 

16eoa^ 


E 


* ^ 

{miai)^(niPi) \J ^g,ni'^g,mi 

r* (z) 

^ nipi,miai\^J i .(j) \2 A{i) 

-I cr,- I ^rriA c 




E 


2r* 

nipi^miCTipiQ. 


CTiPiPi 


\A{'i) 124 ( 4 ) 

\‘^PiQi\ -^rriiGi 


{PiQi) (miCTi) '^g,Pi\l'^g,mi'^g,ni 

Pi > mi 

^ 2r®'? (z) 

E \ " I ^(j) |2 4(i) 

/ ^ • i~~_ ^ I PjQj I mi<Ti 


j = 1 PjQj '^g,Pj y '^g,mi'^g,ni 
3 / i 


E 

PjQkmi 


(z) 

giA/3„-pj. 2 riiPi ,pj Bj qkTkmiai\^J 


X Aid 

Pj Qj Qk Tk mi ai 


j^kjtl=Ai 


>, i = l,..., 4 , ( 29 ) 


where Qjl^pj ^kl^qk ~i~ ciPrni Pijdni 1® the 

wavevector mismatch. 

The coefficients and represent the 

wavevector shift of the optical mode [rii, pi) and the 
linear coupling constant between modes (ni,pi) and 
{mi,ai), induced by the linear perturbations, respec¬ 
tively, and r),.p._„.^.p.p. describe SPM and 

XPM-induced coupling between modes with the same 
frequency, Wi, respectively, PJf.p. mio-ip^gj represents the 
XPM-induced coupling between modes with frequencies 
w* and Wj, and is related to the FWM 

interaction among the pulses. All these nonlinear coef¬ 
ficients have the meaning of ^-dependent effective cubic 
susceptibilities. The linear and nonlinear coefficients in 
Eqs. (l29l) are given by the following relations: 




ojia 

AWL 


[(5efc(r) -b (5eioss(r)]|e„,pAwi)l «‘S', 


riipi ^micr. 


(z) 


^^WLWL, 


(30a) 


X / [(5efc(r) -b (5eioss(r)]e* .(wi) • (30b) 

Js 


riipi ,rnicr- 




2 4 i(cri/3mi-Pi/3ni)2 
{z) = ° p _ _ 


(30c) 


■ riipi imiGiPi Qi 


{z) = 




IT* A/IF* IF* '’s 

^ Pi \J 4ni rii 




■ (^i, -^i,^i)-ep,gi (w*)ep {uJi)emi,Ti {oJi)dS, 


(30d) 




riiPi ,miGiPj Qj 


{z) = 


£2Q4gi(cri/3m.-pi/3p.)z 


WljJWLWL 


f <ipii^^) 

JS 


X^ iOj,Ui).ep.g^[LOj)e [u3j)erniai(!-Ji)dS, 




niPi,PjgjqkTkmiai 


{z) = 


2„4 


eno 




(30e) 

/ KiipSAi) 

Js 


■X^Ag iu}j)e* {u3k)emiai (w/) dS. 


(30f) 


While Eqs. (l29l) seem complicated, in cases of prac¬ 
tical interest they can be considerably simplified. To 
be more specific, these equations describe a multitude 
of optical effects pertaining to both linear and nonlinear 
gratings, including linear coupling between modes with 
the same frequency, nonlinear coupling between modes 
with the same frequency, due to SPM and XPM effects, 
XPM-induced coupling between modes with different fre¬ 
quency, and FWM interactions. In most experimental 
set-ups, however, not all these linear and nonlinear ef¬ 
fects occur simultaneously as in a generic case not all of 
them lead to efficient pulse interactions. 

These ideas becomes clear if one inspects the exponen¬ 
tial factors in Eqs. (j30bl) - (l30eF Thus, they vary over a 
characteristic length comparable to the lattice constant 
of the PhC, namely much more rapidly as compared to 
the spatial variation rate of the pulse envelopes. As 
a result, except for the mode {ni,pi), these linear and 
nonlinear coefficients cancel. There are, however, par¬ 
ticular cases when some of these interactions are phase- 
matched and consequently are resonantly enhanced. To 
be more specific, the integrals in Eqs. (I30bl) - (l30el) are 
periodic functions of z, with period a, so that it is pos¬ 
sible that a Fourier component of these integrals phase- 
matches a specific linear or nonlinear interaction between 
modes (e.g., the linear coupling between two modes with 
the same frequency and SPM- or XPM-induced nonlin¬ 
ear coupling between modes). In this study, we do not 
consider such accidental phase-matching of mode inter¬ 
actions. With this in mind, we discard all terms in 
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Eqs. (|2^ that average to zero to obtain the hnal form 
of the coupled-mode equations for the pulse envelopes: 


{z) ] 51^ iz)P 2 ,n, 52 


9z vl,ni dt 


-:- _•, 


Tlipi 


""g.rii 


2 df^ 

2nvl^^ ^ 


*C<,(2) 




(PiQi) / {rtiPi) 

Pi > rii 




j = 1 Pi Bj 
3 ¥= i 


I \ (z) 

' lniPi,PigiqkTkmi(Ti\^J 


Pjqkrni 

Bj Tfc <^l 

aU) aW* aO) 

Pi Bj BkTk rnitTi 




= 0, i = l,...,4, (31) 


where the new parameters introduced in this equation 
are defined as: 

= K,pA^i)\^dS, (32a) 


j / \ 3w^a2 1 

7n.p.l^j - 


^ / <p.(w*) 


• (uJi, -uJi, uJi):e„.p. (uJi)e* (ujije^.p, (uJi)dS. 


(32b) 


_ 3uJieoa^ 

'yniPi,pipA^) - ^Q i i 

^^g,ni'^g,Pi ^ Snl 




• (Wj, -LOi, 0Ji):ep.g. {oji)e* (a;i)e„,p, {uJi)dS, 


(32c) 


friiPi.pjQj 


{z) = 


iujieoa^ 


16u*_„.u^,p, is„i 






(32d) 


iniPi.pjQjQkTkmiai 


{z) = 


Suiieoo^ 




— —:—-—:—r/ <ip.(wi) • x^^Hwi,-wfc,wi) 

(w{„Wg\w}^^w^X 


■^PjBj (^t)®gtTi, {^k)Gmi(Ti {uJl)dS. 


(32e) 


is much larger than the lattice constant, a. Importantly, 
the linear and nonlinear effects in Eq. appear as be¬ 
ing inverse proportional to the Vg and Vg, respectively. 
In other words, one does not need to rely on any phe¬ 
nomenological considerations to describe slow-light ef¬ 
fects, as they are naturally captured by our model. 


D. Carriers dynamics 


The last step in our derivation of the theoretical model 
describing FWM in Si-PhCWGs is to determine the influ¬ 
ence of photogenerated FCs on pulse dynamics. To this 
end, we first find the rate at which electron-hole pairs 
are generated optically, via degenerate and nondegener¬ 
ate TPA, and as a result of FWM. More specifically, we 
first multiply Eqs. (EH), after all linear terms have been 
discarded, by An^p ., then multiply the complex conjugate 
of Eqs. (12911 by An}pi , and sum the results over all car¬ 
rier frequencies and modes. The outcome of these simple 
manipulations can be cast as: 


E E E E 


dz 


i—l riipi 


OJi 


f=l riiPi 


X Jrri'^ 


r* 

riiPi ,mia. 




2r 


riiPi^miCTiPiQ 


lAii) 12^(7 Aii)* 

\-^mi(7i\ -^miGi-^^niPi 


E 

rriiGi 

(PiBi) ^ {mia-i) '*^g,Pi sj g,mi^g,ni 

Pi > rrii 

^ 2rb' (z) 

\ ^ niPi,miaiPjgj\ E x (7) 124(1) 

Z_^ 1 l~. -^- ' PjBj\ ^miO-i^riiPi 


.\ A (. i ) 124(0 4(0* 

I Vi Pi ' rriiGi riiPi 


j=l PjQj '^g^Pj \l'^g,mi^g,ni 

j * 


(z) 

^iA0r^,VJaum, z niPi,pjgjqkTkmiGi \^) 


PjQkmi 
Qj Tk (71 


<P3^iqkK,myg,n, 


w A ( k )* a { 1 ) 4(0* 

Pj Bj Qk Tk rni ai p. 




(33) 


The sum in the l.h.s. of this equation represents 
the rate at which optical power is transferred to FCs. 
This power is absorbed by carriers generated in the sili¬ 
con slab, in the infinitesimal volume dV{z) = Ani{z)dz, 
where A-a\[z) is an effective area. This area is defined 
in terms of the Poynting vector of the field propagating 
inside the silicon slab. 


In these equations, Sni{z) is the transverse surface of 
the region filled with nonlinear material. Note that the 
exponential factor in the term describing the FWM does 
not average to zero because the FWM interaction is as¬ 
sumed to be nearly phase-matched and therefore the ex¬ 
ponential factor varies over a characteristic length that 


[/^_^J(E(r,t)xH(r,t))t|d^]" 
4J(E(r,t)xH(r,t)),pd5 ' 


(34) 


In this equation, {f)t means the time average of /. Using 
Eq. (001) . and taking into account the fact that An}pi 
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varies in time much slower than e one can express 
Eq. (IMl) in the following form; 


^ 111 ( 2 ;) — 


LJSnl 


EE 

i=l riiPi 


|4 


-lHe(e, 


riipi X 


dS 


4 14(0 |2 


' Snl 


i=l riiPi 


dS 

(35) 


are the same, wi = 0)2 = Wp, whereas the two gener¬ 
ated pulses, the signal and the idler, have frequencies 
0)3 = bjs and W 4 = Wj, respectively. Moreover, we assume 
that all modes are forward-propagating modes and that 
at each carrier frequency there is only one guided mode 
in which the optical pulses that enter in the FWM pro¬ 
cess can propagate - others, should they exist, would not 
be phase-matched - so that we set Ni = 1, i = 1,..., 4. 
Under these circumstances, Eqs. m and Eq. (l36l) can 
be simplified to: 


In spite of the fact that it might seem difficult to use 
this formula to calculate the effective area, we will show 
in the next section that in cases of practical interest it can 
be simplified considerably. We also stress that Eq. (l35ll 
gives the effective transverse area of the region in which 
FCs are generated, so that it should not be confused with 
the modal effective area. In fact, since in the FWM pro¬ 
cess there are several co-propagating beams, a single ef¬ 
fective modal area is not well defined. 

The energy transferred to FCs when an electron-hole 
pair is generated via absorption of two photons with fre¬ 
quencies uii and uij is equal to -I- u)j). Using this 
result and neglecting again all terms in Eq. (I55|) that av¬ 
erage to zero, it can be easily shown that the carriers 
dynamics are governed by the following rate equation; 


dN N 1 


dt Tc hAn\{z) 


EE 

i—\ mpi 






(PiSi) / (riiPi) 
Pi > rii 




j = 1 PjQj 
j / i 


<-Vi -f UJj 


E 

PjQkmi 

Qj Tkcri 


Jm 


Z ^'yniPi,pjQjqkThmiai 

dJi -f Wfc 


V 40 4('=)* 4(0 4(i)* 

PjSj <lk Tfc ^mi a I Pi 




(36) 


where Tc ~ 500 ps [ 6 ^ is the FC recombination time in 
Si-PhCWGs and (' {(") means the real (imaginary) part 
of the complex number, 


IV. DEGENERATE FOUR-WAVE MIXING 

The system of coupled nonlinear partial differential 
equations, Eqs. (1^ and Eq. (IMl) . fully describes the 
FWM of optical pulses and FCs dynamics and repre¬ 
sents the main result derived in this study. In practi¬ 
cal experimental set-ups, however, the most used pulse 
conhguration is that of degenerate FWM. In this partic¬ 
ular case, the optical frequencies of the two pump pulses 


dAp Sp{z)dAp'\ 6 p{z)/d 2 ,pd‘^Ap 

O Q-f-2 


Vg^p dt 
ujpSnfcKpjz) 


5nfcKp{z) ^ , icKp{z)^^ , ^ 3 

■^p (o^fc H” 

nvg^p 2nvg^p 

'piz)\Aj '2 ■ " 


Ip 

ij It' 

i-pl 'T‘2')ps{z)\As\'^ + 2^pi{z)\Ai\^^ Ap 

‘^'5"7p,,(z)A«A,A; = 0, (37a) 

( 5 ^( 2 ) 6s{z)d2,s d'^As 

+ - 2 dt^ 


Vg,s dt 
fcksjz) 


, lCKs(z) , , . 

'As H- - -(ckfc -f (XirdiAs 

^g,S 2flVg_^s 

)|AsP H- 2^sp{z)\-^p\‘^ + 


’-^sp^{z)AlA* = 0, 

k{z) dAj ' 


'-'5 

Si{z)d2,i d'^A, 

2 dfi 


(37b) 


Vg,i dt 

^i{z) ^ , icR^{z) ^ ^ 

- Ai + — -(afc -f Q!in)21i 

,2 

+ ‘2')ip{z)\Ap\^ -f 27 is( 2 :)|Asp] Ai 

lips{z)AlA: = 0, (37c) 


^ U—P-S 


7"(^) 


■\Ap 


z) 


- ■^3m[2Jps^{z) 

UJp 


^ . CJ; 
p—p,s,i ^ 

'fsp^{z)+-f^ps{z)]AlA:A:e-^^^^]}, 

(38) 

- 2/3p. The coefficients of the linear 
in Eqs. (I37l) and Eq. (l38t are: 


u:p)\^dS, 


(39a) 


p)e*p{(^p)ep{ujp)dS, (39b) 

JSni 




el{ujy)ep{ujp)dS, 


(39c) 
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Ipsxiz) = 


SwpEoa^ 


l6Vg^p{Vg^sVg,i)^ Wp{WsWi)^ 


^ ^2 )-^S (^S )®p(^p)^2 (^2 

JSnl 

(39d) 

^ _ SuJseoa'^ _ 1 

16 ^;p,p(^;g,«^;g, 2 )^ Wp{Wsm)^ 

X / e*(a;s) • X^^n‘*^P!-‘^o‘^p):ep(a;p)e*(a;i)ep(wp)d5', 

(39e) 


'S„i 


7ips(^) = 


SuJieoa^ 


l6vg^pivg,sVg,^)^ WpiWsm)^ 

X / e*(a;i) • x^^Hwp,-W2<,Wp);ep(a;p)e*(a;s)ep(a;p)dS', 

(391) 


' S„i 


where ^ and v ^ fi take one of the values p, s, and 
i and the frequency degeneracy at the pump frequency 
has been taken into account. Note that, as expected, 
when the nonlinear coefficients q’s are real quantities, 
namely when nonlinear optical absorption effects can be 
neglected, the optical pumping term in Eq. (1381) vanishes. 
Moreover, since in experiments usually Pp » Ps,Pi, the 
effective area given by Eq. (l35|) can be reduced to the 


a) 


(0 


(d 



- b) - 

(jD 





CO CO 

n V 



0 ) 


ca 


I HV 


c) 



d) 



"^ips’ "^spi 


FIG. 4. Energy diagrams representing the nonlinear opti¬ 
cal processes included in Eqs. ([57)1 . (a) SPM and degenerate 
TPA corresponding to 7 g and jg, respectively, (b) XPM and 
XAM corresponding to 'y'g^ and respectively. Two possi¬ 
ble ways of energy transfer that can occur during a degenerate 
FWM process: (c) two pump photons generate a signal and 
an idler photon, a process described by 7 psi; (d) the reverse 
process, described by 7 ips and ')spi, in which a signal and an 
idler photon generate two pump photons. 


following simplified form: 


^ 111 ( 2 ) 


0 

f 

Snl 

IHe [ep(wp) X h*(a;p)] 

1,25)^ 

J 

f 

Snl 

[gp(u-'p) X 

\^dS 


(40) 


The types of nonlinear interactions incorporated in our 
theoretical model described by Eqs. (E3 are summa¬ 
rized in Eig. |4] via the energy diagrams defined by the 
frequencies of the specific pairs of interacting photons. 
Thus, as per Fig. SKa), the terms proportional to the 7 ^ 
and 7 " coefficients describe SPM and degenerate TPA 
effects, respectively, whereas Fig. SKb) illustrates XPM 
and XAM (also called nondegenerate TPA) interactions 
whose strength is proportional to 7 ^^ and 7 "^^, respec¬ 
tively. Finally, there are two distinct types of FWM pro¬ 
cesses, represented in Fig.|4Kc) and Fig.|4Kd). In the first 
case two pump photons combine and generate a pair of 
photons, one at the signal frequency and the other one at 
the idler, a process described by the term proportional 
to 7 psi. The reverse process, represented by the jips and 
7 spi terms, corresponds to the case in which a signal and 
an idler photon combine to generate a pair of photons at 
the pump frequency. 

As Eqs. (15^ show, the linear and nonlinear optical 
coefficients of the waveguide depend on the index of re¬ 
fraction of silicon, both explicitly and implicitly via the 
optical modes of the waveguide. In our calculations the 
implicit modal frequency dispersion is not taken into ac¬ 
count because it cannot be incorporated in the PWE 
method used to compute the modes. On the other hand, 
the explicit material dispersion is accounted for via the 
following Sellmeier equation describing the frequency de¬ 
pendence of the index of refraction of silicon |69j : 

n (^) - e + ^ + ^2 _ ^2 > 


where Ai = 1.1071pm, e = 11.6858, A = 0.939 816pm^, 
and B = 8.104 61 x lO-^. 

The system of coupled equations, Eqs. (EZl) and 
Eq. (1551) . form the basis for our analysis of degenerate 
FWM in silicon PhC waveguides. In our simulations, 
based on numerical integration of this system of equa¬ 
tions using a standard split-step Fourier method com¬ 
bined with a fifth-order Runge-Kuta method for the in¬ 
tegration of the linear, carriers dependent terms, the z- 
dependence of the coefficients in these equations is rig¬ 
orously taken into account. However, one can signifi¬ 
cantly decrease the simulation time by averaging these 
fast-varying coefficients over a lattice constant, as this 
way the integration step for the resulting, averaged sys¬ 
tem can be increased considerably. The derivation of this 
averaged model is presented in the Appendix. 

One of the key differences between our theoretical 
description of FWM processes in Si-PhCWGs and the 
widely used models for FWM in waveguides with uni¬ 
form cross-section, such as optical fibers or silicon pho¬ 
tonic wires, is that the linear and nonlinear waveguide 
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FIG. 5. (a) Dependence of An\ on z, determined for the odd 
(solid line) and even (dashed line) modes for several values of 
the group-index, Ug. (b) Frequency dispersion of A^i calcu¬ 
lated for the two modes, in the spectral domain where they 
are guiding modes. 



z [a] 



FIG. 6. (a) Dependence of k on z, determined for the odd 
(solid line) and even (dashed line) modes for several values of 
the group-index, rig. (b) Frequency dispersion of R calculated 
for the two modes, in the spectral domain where they are 
guiding modes. 


coefficients are periodic functions of the distance along 
the waveguide. In what follows, we discuss this feature of 
the FWM in more detail, starting with the effective area, 
Ani, defined by Eq. (l40l) . The dependence of this area 
on the longitudinal distance, z, is presented in Fig. [Sja), 
where z spans the length of a unit cell. As we have dis¬ 
cussed, a physical characteristic of slow-light modes is 
their increased spatial extent. This property is clearly 
illustrated in Fig. [SKa), which shows that in the case of 
the even and odd modes the effective area increases by 
almost a factor of two when the group index varies from 
14 to 120 and from 8.6 to 65, respectively. This property 
is also illustrated by the frequency dispersion of the ef¬ 
fective area, averaged over a unit cell, as per Fig. [5Kb). 
Thus, it can be seen in this figure that the effective area 
has a maximum at kz ~ 0.3(27r/a) for the even mode and 
at the edge of the Brillouin zone for both modes, namely 
in the regions of slow light indeed. 

The z-dependence of the spatial mode overlap, k, and 
the frequency dispersion of its spatial average over a unit 
cell, it, are plotted in Figs. |HKa) andllKb), respectively. 
These hgures show that the mode overlap varies more 
strongly with z in the case of the even mode, whereas 
in both cases the mode overlap variation increases as the 
group-index, n^, increases. Interestingly enough, the av¬ 
eraged overlap coefficient of the even mode has a max- 



FIG. 7. Dependence of 5 on z, determined for several values 
of Ug. Solid and dashed lines correspond to the odd and even 
mode, respectively. 


imum at kz Ri 0.3(27r/a), i.e. A ~ 1.52pm, which coin¬ 
cides with a minimum of its Vg. Note also that whereas 
k(z) can be larger than unity within the unit cell, its av¬ 
erage, it < 1. This result is expected because k quantifies 
the mode overlap with the slab waveguide. 

In Fig. [7] we present the dependence on z of another 
physical quantity that characterizes the linear optical 
properties of the PhC waveguide, namely its dispersive 
properties. This parameter, (5, quantifies the extent to 
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A [|jm] A [|jm] 


FIG. 8 . (a), (b) Dependence of 7 ^ and 7 ” on respectively, 
determined for kz = 0.35(27r/a). (c), (d) Frequency dispersion 
of spatially averaged values of 7 ^( 2 ) and Y'{z), respectively, 
determined both for the even and odd modes. 


which the 2 -dependent dispersion coefficients differ from 
their averaged values. Similarly to the mode overlap co¬ 
efficient K, S{z) shows a more substantial changes with 
z in the case of the even mode as compared to the odd 
one and an increase of the amplitude of these oscillations 
with the increase of Ug. Moreover, as it has been demon¬ 
strated in the preceding section, the average of 5{z) over 
a unit cell is equal to unity. 

The 2 -dependence of the nonlinear waveguide coeffi¬ 
cient that characterizes the strength of SPM and TPA 
effects and the wavelength dependence of its average over 
a unit cell are plotted in the top and bottom panels of 
Fig. [51 respectively. One relevant result illustrated by 
these plots is that the nonlinear waveguide coefficient 
increases considerably as the GV of the optical mode 
is tuned to the slow-light regime. Indeed, for the case 
presented in Figs. ma) and | 8 jb) the group-index of the 
odd and even modes are Ug = 6.4 and Ug = 11.5, re¬ 
spectively. This phenomenon is better illustrated by the 
wavelength dependence of the spatially averaged values 
of 7 ^( 2 ) and 7 "( 2 ), which are shown in Figs. [SJc) and 
Hd), respectively. Thus, these plots indicate that the 
nonlinear waveguide coefficient increases by more than 
an order of magnitude as the wavelength is tuned from 
the fast-light to the slow-light regime, the nonlinear in¬ 
teractions being enhanced correspondingly. 


V. PHASE-MATCHING CONDITION 
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FIG. 9. (a), (b) Wavelength diagrams defined by Eq. (14211 and 
Eq. (1441) . respectively. In both panels dashed lines correspond 

to UJp — • 


linear pulse interaction can be achieved. In particular, 
efficient FWM is achieved when the interacting pulses 
are phase-matched, namely when the total (linear plus 
nonlinear) wavevector mismatch is equal to zero. In the 
most general case this phase-matching condition depends 
in an intricate way on the peak power of the pump, 
Pp, signal, Ps, and idler. Pi, as well as on the linear 
and nonlinear coefficients of the waveguide [70l| . This 
complicated relation takes a very simple form when one 
considers an experimental set-up most used in practice, 
namely when the pump is much stronger than the signal 
and idler, Pp :s> Ps,Pi- Under these circumstances, the 
phase-matching condition can be expressed as: 

2j'pPp - 2Pp + / 3 , + / 3 , = 0 . ( 42 ) 

In order to determine the corresponding wavelengths 
of the optical pulses, this relation must be used in con¬ 
junction with the energy conservation relation, that is 
2u!p = ujg OJi- 

An alternative phase-matching condition, less accurate 
but easier to use in practice, can be derived by expand¬ 
ing the propagation constants, j3s,i{oj), in Taylor series 
around the pump frequency, ujp\ 


Before we use the theoretical model we developed to 
investigate the properties of FWM in Si-PhCWGs, we 
derive and discuss the conditions in which optimum non- 




(w - Wp)" 


n>0 


duj''- 


( 43 ) 
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Inserting these expressions in Eq. (1421) and neglecting all 
terms beyond the fourth-order, one arrives to the follow¬ 
ing relation: 

+ /32p(Aw)^ -I- —/34p(Au;)^ = 0, (44) 

where Aw = \ujp — Ws| = |wp — Wi|. 

The wavelength diagrams presented in Figs. IHKa) and 
Kb) display the triplets of wavelengths for which the 
phase-matching conditions expressed by Eq. (|^ and 
Eq. (l4^ . respectively, are satisfied. These wavelength di¬ 
agrams were calculated only for the even mode because 
only this mode possesses spectral regions with anoma¬ 
lous dispersion [cf. Fig. EKb)], which is a prerequisite 
condition for phase-matching the FWM. More specifi¬ 
cally, efficient FWM can be achieved if the pump wave¬ 
length ranges from Ap = 1.52 pm to 1.56 pm. Moreover, 
the diagrams in Fig. IH] show that the predictions based 
on Eq. and Eq. (IT41) are in good agreement, espe¬ 
cially when Aw is small. They start to agree less as 
Aw increases because the contribution of the terms dis¬ 
carded when the series expansion of /3s_i(w) is truncated 
increases as Aw increases. 

Figure |9] also suggests that the spectral domain in 
which efficient FWM is achieved depends on the pump 
power, Pp. To be more specific, it can be seen that 
foi' Pp ^ 0.7 W, a spectral gap opens where the phase¬ 
matching condition cannot be satisfied. The spectral 
width of this gap increases when Pp decreases as the the 
waveguide was not designed to possess phase-matched 
modes in the linear regime. Moreover, the diagrams pre¬ 
sented in Fig. |9] show that in the fast-light regime the 
wavelengths defined by the phase-matching condition de¬ 
pend only slightly on Pp, whereas a much stronger depen¬ 
dence is observed when the wavelengths of the signal and 
idler lie in slow-light spectral domains. 

VI. RESULTS AND DISCUSSION 

In this section we illustrate how our theoretical model 
can be used to investigate various phenomena related to 
FWM in Si-PhCWGs. In particular, we will compare the 
pulse interaction in slow- and fast-light regimes, calculate 
the FWM gain, and investigate the influence of various 
waveguide parameters on the FWM process. The choice 
of the values of physical parameters of the co-propagating 
pulses and that of the input pump power has been guided 
by the exact phase-matching condition given by Eq. (HU). 
In all our calculations we assumed that the pulses prop¬ 
agate in the even mode and, unless otherwise specified, 
the following values for the pulse and waveguide param¬ 
eters have been used in all our simulations: the input 
peak pump power, Pp = lO^Ps = 5W, the input pulse 
width, Tp = Ts = 7 ps, and the intrinsic waveguide loss 
coefficient, = 50dBcm“^. 

Let us consider first the evolution of the envelopes of 
the pulses in the time domain, both in the slow- and 



FIG. 10. Pulse evolution in the time domain. Left (right) 
panels correspond to fast-light (slow-light) regimes, the group- 
index of the pulses being: Ug^i — 9.48 (rig^i = 20.3), ng,p = 
8.64 (ng,p = 8.69), ng,s = 10.37 (ng,s = 23.3). 

fast-light regimes, as illustrated in Fig. [TOl The triplet 
of wavelengths for which the phase-matching condition 
is satisfied is Ap = 1554 nm. As = 1536 nm, and Xi = 
1571 nm in the fast-light regime, whereas in the slow-light 
regime the wavelengths are Ap = 1559 nm. As = 1524 nm, 
and Xi = 1597nm. We stress that in both cases the 
pump pulse propagates in the fast-light regime, whereas 
the signal and idler are both generated either in the fast¬ 
er slow-light regime. 

Under these circumstances, one expects that the pump 
evolution in the time domain is similar in the two cases, 
a conclusion validated by the plots shown in Figs. iinic) 
and frol dl. However, the dynamics of the signal and idler 
are strikingly different when they propagate in the slow- 
light or fast-light regimes. There are several reasons that 
account for these differences. First, whereas the FCA 
coefficient, Ofc, has similar values in the two cases, the 
FCA and intrinsic losses are much larger in the slow-light 
regime because the strength of both these effects is in¬ 
verse proportional to Vg. This is reflected in Fig. [TUI as 
a much more rapid decay in the slow-light regime of the 
signal and idler pulses. Second, it can be seen that in the 
slow-light regime the idler pulse grows at a faster rate. 
This is again a manifestation of slow-light effects. In par¬ 
ticular, the nonlinear coefficient jips, which determines 
the FWM gain, is inverse proportional to [see 

Eq. (13911) ]. As a consequence, the FWM gain is strongly 
enhanced when both the signal and idler propagate in 
the slow-light regime. 

The slow-light effects are reflected not only in the char¬ 
acteristics of the time-domain propagation of the pulses 
but they also affect the evolution of the pulse spectra. 
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FIG. 11. From top to bottom, the left (right) panels show 
the evolution of the spectra of the idler, pump, and signal in 
the case of fast-light (slow-light) regimes. The waveguide and 
pulse parameters are the same as in Fig. 1101 


In order to illustrate this idea, we plot in Fig. [TT] the 
z-dependence of the spectra of the pulses. Similarly 
to the time-domain dynamics, the spectra of the pump 
are almost the same in the slow- and fast-light regimes 
as its GV does not differ much between the two cases. 
The most noteworthy differences between the slow- and 
fast-light scenarios can again be observed in the case of 
the idler and signal. Thus, in the slow-light regime the 
idler decays faster due to increased losses and grows and 
broadens more significantly because of enhanced FWM 
gain and FCD effects, respectively. The influence of FCD 
on the spectral features of the pulses can also be seen 
in the case of the signal and, to a smaller extent, the 
pump. More specifically, Eq. (I^a| shows that the index 
of refraction of the waveguide decreases due to the gen¬ 
eration of FCs. This in turn leads to a^hase-shift and, 
consequently, a blue-shift of the pulse [l^. Interestingly 
enough, one can also see in Fig. [TT] that as the frequency 
of the pulses shifts during their propagation new spectral 
peaks are forming at the initial wavelengths for which the 
phase-matching condition was satisfied. 

In order to gain a deeper insight into the influence 
of slow-light effects on the FWM process, we computed 
the 2 -dependence of the pulse energies when the frequen¬ 
cies of the signal and idler were tuned in the slow-light 
regions of the even mode of the waveguide. We con¬ 
sidered two scenarios, namely these energies were calcu¬ 
lated by including FWM terms in Eqs. (1571) and Eq. (1551) 
and, in the other case, by setting them to zero, that is 
Ipsi = Ispi = lips = 0. In the former case, the EWM 




FIG. 12. (a) FWM enhancement factor vs. propagation dis¬ 
tance, determined for different values of the walk-off param¬ 
eter, A = l/vg.s — l/ag.p. (b) Signal energy vs. propagation 
distance, calculated by including FWM terms in Eqs. ([37|) 
and Eq. (|38|l and by setting them to zero. The blue (green) 
curve corresponds to the fast-light (slow-light) regime consid¬ 
ered in Fig. llOr al [Fig. llOl bl]. whereas the remaining triplet of 
phase-matched wavelengths is Ap = 1556 nm. As = 1530 nm, 
and Ai = 1582nm (red). 


terms are responsible for transferring energy from the 
pump pulse to the signal and idler. Therefore, a suit¬ 
able quantity to characterize the efficiency of this en¬ 
ergy transfer is what we call the FWM enhancement 
factor, 77 , which in the case of the signal is defined as 
77s = 10\og[{EsxF — Esx) /Es^m]- Here, Esxf and Esx 
are the signal energies calculated by taking into account, 
in one case, SPM, XPM, and FWM effects, and only 
SPM and XPM terms in the other case (i.e. FWM terms 
are neglected in the latter case), and Eg^m is the input 
energy of the signal. 

The results of these calculations are summarized in 
Fig. [T2| In particular, it can be clearly seen in Fig. fT^ a) 
that the FWM enhancement factor is strongly dependent 
on pulse propagation regime. To be more specific, as the 
signal and idler are shifting in the slow-light regime a 
smaller amount of energy is transferred from the pump 
pulse to the signal. There are two effects whose com¬ 
bined influence leads to this behavior. First, as we dis¬ 
cussed, the pulses experience larger optical losses in the 
slow-light regime and therefore the signal losses energy at 
higher rate. Equally important, as the pulses are tuned 
in the slow-light regime the walk-off parameter. A, de- 
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a) 




FIG. 13. (a) Dependence of loss factor, A, on the propagation 
distance, z, determined for different values of the group-index, 
Ug. (b) Dependence of A on rig, determined for different val¬ 
ues of a. The slopes of the curves corresponding to z = 10 pm 
and 2 = 70 pm are shown in the inset. 


fined as A = l/vg^g — l/t'g.p, increases, meaning that the 
pulses interact for a shorter time and consequently less 
energy is transferred to the signal. These conclusions are 
clearly validated by the results summarized in Fig. ll2l bl. 
where we plot the energy of the signal vs. the propagation 
distance, determined for several values of the walk-off pa¬ 
rameter. In addition, this figure somewhat surprisingly 
suggests that the FWM process is more efficient in the 
fast-light regime, which is again due to the fact that the 
pump and signal overlap over longer time. 

It is well known that in the slow-light regime linear op¬ 
tical effects are enhanced by a factor of c/vg, whereas cu¬ 
bic nonlinear interactions increase by a factor of (c/vgY'. 
For example, FCA and TPA are proportional to and 
respectively. Our theoretical model predicts, how¬ 
ever, that when the mutual interaction between FCs and 
the optical field is taken into account these scaling laws 
can significantly change. This can be understood as fol¬ 
lows: the amount of FCs generated via TPA, A, is pro¬ 
portional to Vg"^ and since FCA is proportional to the 
product Vg^N, it scales with the GV as 

In order to validate this argument we have determined 
the optical total loss experienced by a pulse when it prop¬ 
agates in the presence of TPA and FCA or, in a different 
scenario, when only TPA is present [the latter case is re¬ 
alized by simply setting = 0 in Eqs. (E3]- Moreover, 
to simplify our analysis, we consider the propagation of 
only one pulse by setting all parameters describing XPM 


and FWM interactions to zero. Finally, we also reduced 
the input power to 100 mW in order to avoid strong SPM- 
induced pulse reshaping. Under these conditions, the 
effect of the FCA on the pulse dynamics can be con¬ 
veniently characterized by introducing a loss factor. A, 
defined as A = 101og[(£’T — ETF)/Ein], where Etf and 
Et are the pulse energies in the case when both TPA 
and FCA terms are included in the model and when only 
TPA is present, respectively, and Ein is the input en¬ 
ergy of the pulse. The results of these calculations are 
presented in Fig. [T^ 

The variation of the loss factor. A, with the propaga¬ 
tion distance, determined for several values of the GV 
is presented in Fig. (nia). As one would have expected, 
the loss factor increases with the group-index, Ug , which 
is a reflection of the fact that the FG-induced losses in¬ 
crease with the decrease of the GV. One can observe, 
however, that when the propagation distance is larger 
than about 139a the loss factor begins to decrease when 
the GV decreases. This behavior is a direct manifes¬ 
tation of slow-light effects, namely as the frequency is 
tuned to the slow-light regime the optical losses increase 
significantly irrespective of the fact that only TPA is con¬ 
sidered or both TPA and FGA effects are incorporated 
in the numerical simulations. 

This subtle dependence of FG-induced losses on Vg 
is perhaps better reflected by the plots presented in 
Fig. dSl^b). Thus, for several values of the propagation 
distance, we have determined the variation of A with 
the group-index, Ug. Then, by calculating the slope of 
the function A{ng) represented on a logarithmic scale 
one can determine how FC-induced losses scale with Vg 
[cf. the inset in Fig. [T^ bi]. The results of this analysis 
clearly demonstrate that FC losses are proportional to 
which agrees with the predictions of our qualitative 
evaluation of this dependence. We stress that for large 
Ug (i.e., small Vg) the v~^ dependence no longer holds 
at large propagation distance, chiefly because the pulse 
is strongly reshaped in the slow-light regime due to en¬ 
hanced nonlinear optical effects, and thus its peak power 
is no longer exclusively determined by optical losses. 


VII. CONCLUSION 

In conclusion, we have derived a rigorous theoreti¬ 
cal model, which describes pulsed four-wave-mixing in 
one-dimensional photonic crystal slab waveguides made 
of silicon. Our theoretical model rigorously incorpo¬ 
rate all key linear and nonlinear optical effects affecting 
the optical pulse dynamics, including modal dispersion, 
free-carrier dispersion, free-carrier absorption, self- and 
cross-phase modulation, two-photon absorption, cross¬ 
absorption modulation, and four-wave mixing. In ad¬ 
dition, the mutual interaction between photogenerated 
free-carriers and optical field is incorporated in our the¬ 
oretical analysis in a natural way by imposing the con¬ 
servation the total energy of the optical field and free- 
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carriers. Importantly, our theoretical formalism allows 
one to derive rigorous formulae for the optical coefficients 
characterizing the linear and nonlinear optical properties 
of the photonic crystal waveguides, avoiding thus any 
of the approximations that are commonly used in the in¬ 
vestigation of nonlinear pulse dynamics in semiconductor 
waveguides based on photonic crystals. 

As a practical application of the theoretical results 
developed in this study, we have used our theoretical 
model to investigate the properties of degenerate four- 
wave-mixing of optical pulses propagating in photonic 
crystal waveguides made of silicon, with a special focus 
being on highlighting the differences between the pulse 
dynamics in the slow- and fast light regimes. This anal¬ 
ysis has revealed not only that linear and nonlinear ef¬ 
fects are enhanced in the slow-light regime by a factor 
of Ug and n^, respectively, but also that these scaling 
laws are markedly affected in the presence of free-carriers. 
Moreover, since our study has been performed in a very 
general framework, i.e. generic optical properties of the 
waveguides (multi-mode waveguides) and pulse config¬ 
uration (multi-frequency optical field), our findings can 
also be used to describe many phenomena not considered 
in this work. For example, important nonlinear effects, 
including stimulated and spontaneous Raman scatter¬ 
ing, coherent anti-Stokes Raman scattering, and third- 
harmonic generation, can be included in our model by 
simply adding the proper nonlinear polarizations. 
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APPENDIX: AVERAGED MODEL DESCRIBING 
DEGENERATE FOUR-WAVE-MIXING 

The spatial scale over which the envelope of picosec¬ 
ond pulses varies is much larger than the lattice constant 
of the PhG and therefore for such optical pulses one can 
simplify the system of equations governing the pulse in¬ 
teraction, i.e. Eqs. dSZl) and Eq. (1551) . by taking the av¬ 
erage over one lattice constant. Under these conditions, 
the corresponding system of coupled equations can be 
cast in the following form: 

. / dAp I (9Ap \ /^2,p ^ Ap ujpdfi^cKp . 

92 Vg^p dt ) 2 dt'^ ^ 

+ „ (ofc + Oin)Ap -|- 2e*^^^7psi 
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A A 27 ps|A;,p -|- 27 pi|Aip) Ap = 0, (45a) 
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(46) 

The coefficients of the linear and nonlinear terms in these 
equations are given by the following formulae: 
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FIG. 14. Comparison between pulse evolution as described 
by the full and averaged model is presented in the left and 
right panels, respectively. The group-index of the pulses are 
Ug^i = 20.3, ng,p = 8.69, and ng_s = 23.3 and correspond to 
the slow-light propagation scenario presented in Fig. 1101 The 
bottom panel shows the ^-dependence of the normalized pulse 
amplitude, 4'p(2) = Aij,{zq +z) / Ag,{zo), fi —p,s,i, calculated 
for the unit cell starting at zo = 200a. 


where V^i is the volume occupied by silicon in a unit 
cell of the PhC waveguide and zq an arbitrary distance. 
Finally, An\{z) in Eq. ( |47gP is given by Eq. (l40l) . whereas 
the index vo takes any of the following values: p, s, i, ps, 
si, ip, psi, spi, or zps.whereas the index va takes any of 
the following values: p, s, i, ps, si, ip, psi, spi, and ips. 
Note that in deriving the averaged equations governing 
the pulse and EC dynamics, Eqs. (H51) and Eq. (H5)) . we 
have assumed that the FWM process is nearly phase- 
matched, namely 1/a. In other words, the phase 

of the exponential factors in these equations vary much 
slower with the distance, z, as compared to the variation 
of the dielectric constant of the PhC waveguide. 

A comparison between the predictions of the full and 
averaged models is illustrated in Fig. [TH Thus, we have 
considered the slow-light pulse dynamics presented in 
Fig. [To] and determined the pulse evolution using both 
the full and averaged models. As it can be seen, both 
models predict a similar pulse dynamics for the entire 
propagation length, z = 1000a. This result is expected 
as the envelope of picosecond pulses, as are those chosen 
in our simulations, spans a large number of unit cells and 
therefore the pulse amplitude is only slightly affected by 
the local inhomogeneity of the index of refraction. 

The fast variation with z of the pulse envelope is shown 
in Fig. fTTl' g'l. where we plot the z-dependence of the nor¬ 
malized pulse amplitude, 4'^(z) = A^(zo + z)/A^(zo), 
pL = p,s,i, calculated for the unit cell starting at zq = 
200a. It can be seen in this figure that the pulse en¬ 
velope varies at a spatial scale commensurable with the 
lattice constant yet the amplitude of these variations is 
much smaller than the pulse peak amplitude. The mag¬ 
nitude of these variations, however, would comparatively 
become more significant should the pulse duration would 
be brought to the femtosecond range. 
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